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Abstract 

In this paperrthe problem of random vibration of ge- 
ometrically nonlinear MDOF structures is considered. 
The solutions obtained by application of two different 
versions of a stochastic linearization method are com- 
pared with exact (F-P-K) solutions. 

The formulation of a relatively new version of the 
stochastic linearization method (energy-based version) 
is generalized to the MDOF system case. AlsoFa new 
method for determination of nonlinear stiffness coef- 
ficients for MDOF structures is demonstrated. This 
method in combination with the equivalent lineariza- 
tion technique is implemented in a new computer pro- 
gram. 

Results in terms of root-mean-square (RMS) dis- 
placements obtained by using the new program and an 
existing in-house code are compared for two examples 
of beam-like structures. 

Introduction 

Resurgent interest in high speed fight vehicles and 
the daily operation of the aging commercial and mil- 
itary aircraft fleets necessitate the further develop- 
ment of sonic fatigue technology to understand the fa- 
tigue mechanisms and to estimate the service life of 
aerospace structures subjected to intense acoustic and 
thermal loads. Efforts to extend the performance and 
flight envelope of high speed aerospace vehicles have 
resulted in structures which may behave in a geomet- 
rically nonlinear fashion to the imposed loads. Such 
behavior can have a significant effect on fatigue life. 
Further improvements in vehicle performance and sys- 
tem design are hampered by the limited understand- 
ing of the physical nature of geometrically nonlinear 
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structural response. Conventional (linear) prediction 
techniques can lead to grossly conservative designs and 
provide little understanding of the nonlinear behavior. 
A large body of work exists on the prediction of ge- 
ometrically nonlinear dynamic response of structures. 
All methods currently in use are typically limited by 
their range of applicability or excessive computational 
expense. 

Methods currently in use to predict geometrically 
nonlinear dynamic structural response include per- 
turbationrFokker-Plank-Kolmogorov (F-P-K)rMonte 
Carlo simulation and stochastic linearization tech- 
niques. Perturbation techniques are limited to weak ge- 
ometric nonlinearities. The F-P-K approach 1,2 yields 
exact solutionsr but can only be applied to simple 
mechanical systems. Monte Carlo simulation is the 
most general methodrbut computational expense lim- 
its its applicability to rather simple structures. Finallyr 
stochastic linearization methods (e.g. equivalent lin- 
earizationFsee 2-6 ) have seen the most broad applica- 
tion for prediction of geometrically nonlinear dynamic 
response because of their ability to accurately capture 
the response statistics over a wide range of response 
levels while maintaining a relatively light computa- 
tional burden. 

Implementations of stochastic linearization have 
been limited to special purpose computer codes un- 
til recently when the method of equivalent lineariza- 
tion was introduced into MSC/NASTRAN as a Direct 
Matrix Abstraction Program (DMAP) Alter 7 . In this 
study an alternative approach to the solution of non- 
linear vibration problems is developed and an inde- 
pendent in-house code based on this approach is im- 
plemented. 

Equivalent Linearization Techniques 

Two versions of the equivalent linearization tech- 
nique are considered. One is based on minimization of 
the error in the force- vectorr and the other minimizes 
the error in potential energy. 

Force Error Minimization Version 

Consider a MDOFFviscously damped linear system. 
The equations of motion governing such a system can 
be written in the form 

MX + CX + KX = F 
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where M is the mass matrixr C is the damping ma- 
trixr K is the stiffness matrixr X is the displacement 
response vector and F is the force excitation vector. 

For geometrically nonlinear problems of deforma- 
tionr e.g. large deflection flexural vibration of thin 
plate structuresr the governing equation(s) of motion 
will include a nonlinear force term r(X)Ti.e. 

MX + CX + KX + r(X) = F (1) 

where the vector function r(X) generally includes 2nd 
and 3rd order terms in X. There exist mathematical 
difficulties in the derivation of a general solution to 
equation (1) for the case of random excitation. An 
approximate solution can be achieved by formation of 
an equivalent linear system: 


MX + CX + (K + K e )X = F (2) 

where K e is the equivalent linear stiffness matrix. 

The method of equivalent linearization seeks to mini- 
mize the difference between the nonlinear force and the 
product of the equivalent linear stiffness and displace- 
ment response vector. The equivalent linear stiffness 
satisfying this requirement can be determined from the 
following condition: 


error = E [(r(X) - I< e X) T (T(X) - K e X)] — >■ min 


where E[..i\ represents the expectation operator. The 
latter equation will be satisfied if 


d {err or) 
dKeii 


= 0 


i,j = 1,2, ...,N 


In this paperFconsideration is limited to the case of 
Gaussianr zero-mean excitation and response to sim- 
plify the solution. Omitting intermediate derivationsr 
the final form for the equivalent linear stiffness ma- 
trix becomes (see for example Roberts et al. 3 rAtalik 
et al. 4 ): 


A-. = (3) 

Potential Energy Error Minimization Version 

Elishakoff et al. 5,6 proposed another stochastic lin- 
earization approachr based on potential energy error 
minimization and numerical results were demonstrated 
for the case of SDOF systems. In this paperrthat ap- 
proach is generalized for the case of MDOF systems. 
One can begin with an expression for the error in po- 
tential energy e 

e = E[{U{X)- l -X T K e Xf] 

where U (A) is the potential energy of the original (non- 
linear system) and K e is the stiffness matrix of the 
equivalent linear system. 


A condition of minimized error e —>■ min requires the 
following 

E[ W~ iiu{x) ~ \ xTReX ^ = 0 *'>•* = 1 > N 

where K e ij are elements of matrix K e . Omiting inter- 
mediate derivationsFone obtains the following system 
of N 2 linear equations with respect to unknown ele- 
ments of matrix K e : 

N 

^ K e ijE[xiXjX k xi\ = 2E[x k xiU{X)\ (4) 

i,j = 1 


where h,l = 1, N . 

For examplerequation (4) would have the following 
form for a two-degree-of-freedom system 


E[x f\ E[x\x 2 ] E[x\x 2 \ E[x\x^\ 


Kell 

E[x\x 2 \ £'[* 1 * 2 ] E[x'fx'{] E[xixf\ 


Ke 12 

E[x\x 2 ] E[x\x^\ E[x\x^\ Af*!* 3 ] 


K e '21 

E[x\x 2 ] Af*!* 3 ] Af*!* 3 ] -^[* 2 ] 


Kg 22 


= 2 


E[x\U{X)} 

E[x 1 x 2 U(X)\ 

E[x 1 X2U{X)\ 

E[x\U{X)\ 


(5) 


Note that the 2nd and 3rd rows are identicairthus an 
additional equation is required to solve this system. 
The additional equation(s) can be provided by the im- 
position of a condition of symmetry of the matrix K e : 


Keij — K e ji (6) 

The matrix of the system in equation (5) involves 
4th order moments of displacements and the right-hand 
side (assuming that the potential energy is a function 
of the 2ndr3rd and 4th order terms) involves moments 
of 4thr 5th and 6th order. Using the Gaussian dis- 
tributedr zero-mean response assumption means that 
the odd order moments are zero and the higher even 
order statistical moments can be expressed in terms of 
the 2nd order momentsre.gT 


E[xiXjX k xi] = E[x i x j ]E[x k x l \ + E[xiX k \E[xjXi\- \- 


E[xiXi\E[xjX k \ 

Therefore the matrix and right-hand side of (5) can 
be determined solely by the response covariance ma- 
trix. So the equivalent stiffness matrix at each iteration 
is determined through the use of response covariance 
terms from the previous iteration by solving equations 
(5). 
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Iterative scheme of equivalent linearization 

Having defined the equivalent linear stiffness matrix 
through either the force error or potential energy er- 
ror minimization techniquesFone can proceed with the 
solution of the equivalent linear system. Assuming sta- 
tionary excitationra stationary response is sought pre- 
cluding the need for initial conditions. As the equiva- 
lent stiffness matrix K e is a function of the unknown 
displacement response vectorrthe solution to the sys- 
tem of equations of motion takes an iterative formri.e. 

MX n+1 + C'X n+ i + (K + K en )X n+1 = F (7) 

where new displacement response estimates are calcu- 
lated from a system based upon the previous estimate 
and iterations are continued until a convergence crite- 
rion is satisfied. 

The solution to the equivalent system in equation 
(7) for each iteration can be obtained in the frequency 
domain using the well known relation between the spec- 
tral density matrices for a linear system: 

5a;(w)„ + i = H n (uj)S f (uj)H T n (uj) 

where the over-bar indicates the complex-conjugaterS'/ 
is the spectral density matrix of the random excitation 
and the frequency response matrix is given by 

H n { ui) = [— ui M + iuiC + K + K en \ 

The zero-time-lag covariance matrix components 
participating in the matrix K en are calculated from 
the response spectral density matrix using the Wiener- 
Khinchine formula 

/ OO 

[Sxij (^)]n duj 

- oo 

An implementation of the equivalent linearization 
approach outlined above was recently implemented in 
a special purpose in-house code to provide a tool for 
expedient study. 

Comparison with F-P-K Solutions 

The two equivalent linearization methods presented 
above will be compared with F-P-K solutions for SDOF 
and 2DOF systems. 

SDOF system 

Consider a SDOF system (Duffing oscillator): 

q(t) + 2 & 0 q(t) + ojlq(t) + /3q 3 (t) = f(t) (8) 

where q is a nondimensional coordinate/displacement. 
The addition to the potential energy originated from 
the nonlinear term is characterized by 

U(q) = \pq A ( 9 ) 


For this caseFsolution of the system (5) (energy-based 
technique) provides the following equivalent stiffness 

k e = 2.5 flE[q 2 ] 

and equation (3) (force-base technique) yields 
k e = 3 f3E[q 2 ] 

Comparison of response variances for this system 
versus the nonlinearity parameter is illustrated 
in Figure 1. A white noise excitation was taken as 
the input spectral density functionr i.e. Sf(ui) was 
constant and equal to l.Oe + 05. The results corre- 
spond to an oscillator with a natural frequency of 
57.4 El z (w 2 = 1.301e-|-05 s -2 ) and damping coefficient 
£ = 0.005. The three curves in Figure 1 correspond 
to the F-P-K solutionr force error minimization and 
energy error minimization versions. 

Comparison of response variances for this oscillator 
versus the spectral density function value Sj is illus- 
trated in Figure 2T where the nonlinearity parameter 
was fixed and equal to 10. One can see that the 

“o ... 

energy error minimization version results are closest to 
the exact (F-P-K) solution results. 

2DOF system 

As a next exampler consider the model of a 2DOF 
system in Figure 3. The equations of motion for this 
model have the form 

mqi + ciqi+kiqi+k 2 (qi-q2)+/3iqf+/32(qi-q2) 3 = fl(t) 

m 2 q 2 + c 2 q 2 + k 2 (q 2 ~ ?l) + P 2 {q 2 ~ qif = f' 2 (t) (10) 

The potential energy contribution from the nonlinear 
part is 

U (qi, q 2 ) = + ^/3 2 (?2 - ?i) 4 

White noise was taken again as the input exci- 
tation with the spectral density matrix components: 
S/u(w) = Sf 22 (u>) = IT and S/i 2 (w) = S/ 2 i(w) = 0. 
The rest of the parameters of this model were as fol- 
lows mi = m 2 = IT ki = k 2 = lTci = c 2 = 0.1 and 
/3i = f3 2 = /3. A comparison of response variances ver- 
sus the nonlinearity parameter /? for this 2DOF system 
(10) is shown in Figure 4. Again the energy-based ver- 
sion results are closer to the exact F-P-K solutionrthan 
the force-based version results. 

Note that in the case of general MDOF nonlinear 
systemsr the determination of the expression for the 
potential energy can be complicated. This problem 
will be addressed in the section below. 
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Determination of Nonlinear Stiffness 
Coefficients 

So faiTexamples were considered where the nonlin- 
ear stiffness coefficients were prescribed. In a general 
case of a MDOF systemrthese coefficients have to be 
determined. One method of determining the nonlin- 
ear stiffness coefficients is through the use of a finite 
element approach. Existing finite element commercial 
programs are unable to provide these nonlinear stiff- 
ness coefficients directly. It is desirable to achieve a so- 
lution within a commercial finite element code to take 
advantage of the comprehensive element libraryr etc. 
necessary for modeling complex structures. This sec- 
tion describes a method of determining the nonlinear 
stiffness coefficients through the use of the nonlinear 
static solution capability that exists in many commer- 
cial finite element codes. 

For MDOF structuresr it is expedient to seek a so- 
lution in modal coordinate space 

X = <S>q (11) 

where $ is generally a subset (L < N) of the linear 
eigenvectors (normal modes). Such a representation 
allows the size of the problem to be significantly re- 
duced without a noticeable loss of accuracy in many 
cases. 

One can obtain the following set of differential equa- 
tions in terms of modal coordinates qi (i = 1, L): 

L 

Qi(t) + ^2 c f?9j(t) + + 7i(9i, 92, •••, q£) = fi(t ) 

j = i 

. (12) 

where the nonlinear terms will be represented in the 
following form 

L L 

7;(9i,92, ...,9l) = X] a )kWk + X] & i«9j'9fc9i 

j,k=j j,k=j,l=k 

. ^ 

where the first index j takes values H2E..ILrthe index 
k takes values from j (the current first index value) 
and up j + II}" + 2 ... to L and the third index l takes 
values from k (the current second index value) and up 
k + lTA; + 2 to L. 

The analytical form of the nonlinear terms facilitates 
the solution of equations (12) when the forces and dis- 
placements are random functions of time. 

A procedure for determination of the coefficients a l - k 
and b) kl is described briefly. This procedure requires 
the application of a finite element program with a non- 
linear static solution capability. In this studyr the 
MSC/PATRAN and MSC/NASTRAN programs 9 - 10 
are utilized. 

The suggested technique is based on the restoration 
of nodal applied forces from enforced nodal displace- 
ments prescribed to the whole structure in a static so- 
lution (linear and nonlinear). Namelyrby prescribing 


the physical nodal displacements (vector X c ) to the 
structureFone can restore the nodal forces Ft and the 
corresponding nonlinear contribution F c : 

F c = T(X e ) = F t - I<X c (14) 

The displacements X c can be prescribed by creating a 
displacement constraint set for the model in PATRANT 
then the nodal applied forces Ft will arise as single- 
point-constraint forces in a NASTRAN nonlinear static 
solution. 

To illustrate the techniquer one can begin with the 
prescription of displacements for the whole structure 
in the following form 

X c = Mi ( 15 ) 

The nodal force vectors Ft (nonlinear static solution) 
and KX C (linear static solution) are provided by NAS- 
TRAN. The nonlinear term F c can then be evalu- 
ated by equation (14). The vector of modal forces 
F c = $ T F c is calculated and it is represented as 

f c = <s> t f c = $ T r(x c ) = $ T r(Mi) = 

Ki]9i9i + Kn]9i9i9i (16) 

where the sought stiffness coefficients [aii]r[& ! m ] are 

column-vectors L x 1 (i = 1 ,L). Note that all other 
nonlinear terms in (16) do not appear since qj = 0 for 

i ± 1- 

Prescribing a displacement field with opposite sign 
X c = —<f>iqi results in a modal force vector (denoted 
by T-c): 

F. c = $ t T_ c = $ T r(x c ) = 

® T T(~Mi) = Ki]9i9i - [&m]9i9i9i (17) 

where the quadratic (even) term will be the same as in 
(16) and the cubic (odd) term takes on a sign change. 

Note that in the system of equations (16) and (17)T 
the value of q\ is given. The coefficients [a* 1:L ] r[6® 11:L ] 
(i = 1, L ) can be determined from this system of 2 xi 
linear equations. In an analogous mannerr i.e. pre- 
scribing X c = (fijqjT all other coefficients [a®- •] r [6®- • •] 
can be determined. 

A similar technique can be employed to determine 
coefficients with two or three inequal lower indicesFe.gT 
M&iiJF [&i 22 ] or [^ 123 ] - Note that coefficients of 
the latter type appear only if the number of retained 
eigenvectors L in (11) is greater than or equal to 3. De- 
termination of coefficients [“ 12 ] ^ 112 ] and [b\ 22 ] will be 
considered as an example. Prescribe the displacement 
field to the model in the following form 

X c = Mi + Ml2 

then the calculated (using NASTRAN) corresponding 
modal force vector F c is represented as follows 

F c = $ t T c = $ T r(c?iigi + M 2 ) = 
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[ a ll]?l?l + [frlllklSlSl + [ a 22]?2?2 + [&222]?2?2?2 + 

K 2 ]?i ?2 + [b\ 12 \qiqiq 2 + [b\ 22 \q 2 q 2 qi (18) 

Prescribing the opposite sign displacement field 

X c = -<f>iqi - </ 2 ?2 
one obtains a second set of equations 

F-c = $ T F_ c = $ T r(-^igi - <f) 2 q 2 ) = 

Kl]?l?l - [&lll]?l?l?l + [a* 22 ]?2?2 - [^222]92?2?2 + 

K 2 ]?i? 2 - [&n 2 ]Si9i92 - [b\ 22 \q2q2qi (19) 

Summing (18) and (19)rone obtains 

F c + F-c = 2[a\ 1 ]qiqi + 2[a l 22 \q 2 q 2 + 2[a\ 2 ]qiq 2 

From this equationr the coefficients [ a\ 2 ] are deter- 
mined (note that the coefficients [ai-jrfa^] were al- 
ready determined above). 

Now we have two sets of L equations (18) and (19)T 
but to determine cubic coefficients [&* 112 ] and [&* 122 ] from 
them is not possible since the system matrix has lin- 
early dependent rows. ThereforeFan additional type of 
displacement Held is required. One can prescribe the 
following type 


x a = <t>iqi - 4>2q2 

Then the modal force vector is equal to 

F a = <P T F a = $ T r((?iigi - <f) 2 q 2 ) = 

Kiki?i + [b\n\qiqiqi + [a 22 ]g 2 g 2 - [b 222 \q 2 q 2 q 2 - 
K 2 ki92 - [b\ 12 \qiqiq 2 + [b\ 22 \q 2 q 2 qi ( 20 ) 

From the system of 2 *L linear equations (18) and (20)T 
the coefficients \b\ 12 \ and [&* 122 ] can be determined. In 
a similar mannerFall coefficients of the type [ 6 ®. -j.] and 
\b\kj\ can be determined. A technique has been devel- 
oped to determine all the coefficients using 

a similar approach as above. 

Solution of modal equations 


Having the modal equations of motion (12) formu- 
latedr solution to these equations can now be under- 
taken through a variety of techniques. For the case 
of random loadingr the application of the equivalent 
stochastic linearization was implemented in this study. 
Within the framework of the force-based techniquerthe 
equivalent stiffness matrix (according to the formula 
(3)) will have the following form 


Ke 


^(71,72, ...,7 l ) 
V d(qi,q 2 , ...,q L ) 


(21) 


Note that the derivatives and expectations in (21) can 
be easily evaluated due to the analytical representation 


of the nonlinear terms in (13). A program producing 
the calculations described above has been developed 
and numerical results will be demonstrated in the next 
section. 

Based upon the expressions derived in equation (13)T 
one can proceed with the determination of potential 
energy U in terms of modal coordinates. It is known 
that elastic force terms (linear + nonlinear) satisfy the 
following 

dU 

kiqi +7i(qi, ? 2 , ■■■,qL)= ^7 i = l,L (22) 

Since all nonlinear coefficients in ji(qi , q 2 , ..., q£) 
have been determinedr the potential energy function 
U(qi, q 2 , ..., (/l) can be derived and it can be used in the 
energy-based stochastic linearization technique. An 
implementation of the energy-based version for appli- 
cation to MDOF systems is considered as future work. 

Numerical Results for MDQF Structures 

It is important to note that the analysis of a vi- 
brating structure in the nonlinear setting is necessary 
only if the comparison of two static solutions (linear 
and nonlinear ones for the highest deformation level) 
shows a noticeable difference in the displacement fields. 
To illustrate thisr one can consider the three beam 
structures shown in Fig.5a-c. It was found that for 
a cantilevered beam model (Fig. 5c)Tthe difference in 
terms of static flexural displacements is negligible. Two 
curves (Fig. 6 ) corresponding to the linear and nonlin- 
ear models are indistinguishableFwhere the flexural dis- 
placement of the tip node versus the applied static base 
acceleration (inertial loading) is plotted. However for a 
beam in Fig. 5b and clamped-clamped beam (Fig. 5a) 
the difference in displacements is quite noticeable (see 
Fig. 6I7)rindicating that the vibration analysis should 
be conducted in the nonlinear setting. 

The numerical results presented in this section cor- 
respond to models of the structures in Figures 5a and 
b. The results obtained with the SEMELRR DMAP 7 
are compared with the new method which employs 
the technique described aboveri.e. the determination 
of nonlinear stiffness coefficients plus the conventional 
(force-based) stochastic linearization technique. 

SEMELRR was implemented in MSC/NASTRAN 
using equivalent linear modal degrees of freedom. This 
requires an eigensolution at each iterationr but af- 
fords the most simple and versatile procedure read- 
ily adapted within the framework of the existing 
MSC/NASTRAN solution sequences. The original im- 
plementation was limited to spatially uniform mechani- 
cal loadsrbut has since been generalized to include spa- 
tial non-uniformity. The solution is also formulated to 
include the effects of static deflection due to mechani- 
cal or thermal loadsFmaterial nonlinearity and follower 
forces. Some work has been done to validate the pre- 
diction capability of SEMELRRFsee Robinson et al. 7 r 
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material 

Young’s 

modulus 

Poisson’s 

ratio 

density 

aluminum 

0.73e+ll 

0.325 

0.2763e+04 

length 
beam a) 

length 
beam b) 

width 

thickness 

0.4572 

0.2286 

0.0254 

0.002261 


Table 1: Parameters of beams a) and b) 


b ln 

^222 

^112 

^122 

0.899e+12 

0.977e+13 

0.191e+13 

0.139e+14 

P 

°111 

^222 

P 

°112 

^122 

0.638e+12 

0.608e+14 

0.139e+14 

0.293e+14 


Table 2: Nonlinear stiffness coefficients for beam a) 

but rigorous validation of all of its features and bounds 
of its applicability have not been fully established. 

The parameters of the models in Fig. 5a and b are 
shown in Table 1 (system of units is Sir[m]r[iV/m 2 ]r 
[kg] etc.)rwhere width and thickness are dimensions of 
the cross-section of the beams. The first two natural 
frequencies (associated with flexural modes in the ex- 
citation plane) for the beam in Fig. 5a are 57.4 Hz and 
310.1 Hz and the first two natural frequencies for the 
beam in Fig. 5b are 35.6 Hz and 220 Hz. In all cases 
the first two (symmetric for clamped-clamped beam) 
flexural mass-normalized modes were chosen to approx- 
imate the motion of beams according to formula (11). 

The nonlinear stiffness coefficients determined with 
application of the procedure described above are sum- 
marized in Tables 2 and 3. The quadratic terms were 
negligibleFso only the 3rd order terms are shown. Since 
the modal coordinates qi,q 2 are nondimensionairthe 
units of these nonlinear coefficients are in [N * m] . 

Note that from (22) would follow that 

djj _ djk_ _ d 2 U 
dq k dqj dq k qj 

and comparing the terms with like powers in qj and 
q k leads to a certain relation between the nonlinear 
coefhcientsrfor examplerfor the cubic coefficients b ] 2 2 
and b\ 12 it is 

b 1 - b 2 

°122 ~ u 112 

and for other typesrit is 

3&222 = ^122 3& 2 n = b\ 12 

It turned out that the computed nonlinear stiffness co- 
efficients (see Tables 2 and 3) are in an excellent agree- 
ment with these relations. 

The results in terms of the RMS displacement of the 
middle and tip nodes (Fig. 5a and b) are shown in 


b ln 

^222 

^112 

b l 

°122 

0.359e+13 

-0.713e+14 

-0.233e+14 

0.722e+14 

P 

°111 

^222 

P 

°112 

^122 

-0.779e+13 

0.247e+15 

0.722e+14 

-0.213e+15 


Table 3: Nonlinear stiffness coefficients for beam b) 

Fig. 8 and 9. A vertical base white noise excitation 
(acceleration a k ) provided inertial loading which was 
spread over a 20-320 Hz range. One can see that nu- 
merical results produced with the SEMELRR code and 
the new method differ by about 20 % for the case of 
clamped-clamped beam. The difference is about 30 % 
for beam b). In each caserthe SEMELRR’s RMS dis- 
placements are less than the RMS displacements from 
the new approach. 

Unfortunatelyrthere are no exact solutions available 
for these structuresr so comparisons are not possible. 
HoweverFrecent experimental measurements (not pre- 
sented here) indicate that the new method predicts 
RMS responses more in agreement with their physical 
counterparts than the SEMELRR solution sequence. 
This will be quantified with further numerical and ex- 
perimental work. 

Summary 

The energy-based version of stochastic linearization 
technique has been extended to MDOF systems and 
the numerical results have shown superior performance 
of this technique in comparison with the conventional 
linearization version. 

A new method for determination of nonlinear stiff- 
ness coefficients has been suggested and applied to sev- 
eral examples of beam-like structures. This method 
has been incorporated into a program which cal- 
culates a steady-state response of a MDOF struc- 
ture to a Gaussian zero-mean excitation. Efforts are 
presently underway to implement this capability into 
MSC/NASTRAN through a DMAP Alter. 

Some difference (about 20-30 % range) has been 
found between the two independent results in terms 
of prediction of nonlinear response. Further numerical 
studies and experimental work will be devoted to this 
problem. 
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Figure 3: Example of a two-degree-of- freedom system 



Figure 2: Variance of displacement vs. spectral density 
of force for the SDOF system with natural frequency 
wo/27t = 57.4 77 :T damping coefficient £ = 0.005 and 
13/loI = 10 
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Variance of displacement 




Figure 6: Displacement as a function of inertial loadr 
Figure 4: Variancies of displacement vs. nonlinearity beams in Fig.5bE 
parameter for the 2D0F system 
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Figure 5: Examples of beam-like structures 


Figure 7: Displacement as a function of inertial loadr 
beam in Fig. 5a 
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Figure 8: RMS displacement. for beam in Fig. 5a. as a. 
function of inertial random loading 



Figure 9: RMS displacement, for beam in Fig. 5b as a. 
function of inertial random loading 
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